Discrete radon transform

This paper describes the discrete Radon transform (DRT) and the exact inversion algorithm for it. Similar to the discrete Fourier transform (DFT), the DRT is defined for periodic vector-sequences and studied as a transform in its own right. Casting the forward transform as a matrix-vector multiplication, the key observation is that the matrix-although very large-has a block-circulant structure. This observation allows construction of fast direct and inverse transforms. Moreover, we show that the DRT can be used to compute various generalizations of the classical Radon transform (RT) and, in particular, the generalization where straight lines are replaced by curves and weight functions are introduced into the integrals along these curves. In fact, we describe not a single transform, but a class of transforms, representatives of which correspond in one way or another to discrete versions of the RT and its generalizations. An interesting observation is that the exact inversion algorithm cannot be obtained directly from Radon’s inversion formula. Given the fact that the RT has no nontrivial one-dimensional analog, exact invertibility makes the DRT a useful tool geared specifically for multidimensional digital signal processing. Exact invertibility of the DRT, flexibility in its definition, and fast computational algorithm affect present applications and open possibilities for new ones. Some of these applications are discussed in the paper. T INTRODUCTION HE discrete Radon transform (DRT) is a discrete version of the classical Radon transform (RT) [l] and some of its generalizations [2]-[4]. The DRT defined and described in this paper is exactly invertible in an efficient manner. Since RT (and, hence, DRT) have no nontrivial analog in the one-dimensional space, exact invertibility makes DRT a useful tool geared specifically for multidimensional digital signal processing. Discrete versions of the classical RT are being used in signal processing and there is an extensive literature devoted to this subject. Procedures which are discrete versions of the RT are known as slant stack [5], tau? transform [6], [7 ] , velocity filtering [8], [9], fan filtration [9], and beam forming [lo]. These procedures were successfully used in various applications such as ground roll removal [8], plane wave decomposition [ 111-[ 141, P-S separation [ 151, interpolation and resampling of data [8], and, also, in procedures such as velocity analysis and beam steering. Some of these applications depend on invertibility of the RT. There are two major classes of algorithms for inversion of discrete versions of the RT: algebraic reconstruction techniques, and backprojection algorithms based upon Manuscript received January 17, 1986; revised July 10, 1986. The author is with Schlumberger-Doll Research, Ridgefield, CT 06877IEEE Log Number 8610882. 4108. various discretizations of Radon’s inversion formula. We note that in applications mentioned above, algebraic reconstruction techniques (such as iterative and row-action methods [ 161-1181 developed for image reconstruction from projections) were not used. Although these procedures can, in principle, provide an exact inversion of the RT, they are not practical in these applications because of the size of the linear system to be solved. As a consequence, all inversion algorithms used, in practice, have been approximate inversions based on Radon’s inversion formula. We show, however, that a discrete version of the RT can be inverted only approximately if its inversion is based on a discretization of Radon’s formula. This inversion formula was discovered by J. Radon [l] who was also the first to define the transform as such. The transform itself and inversion formula were later rediscovered for use in different applications which include astrophysics [19] and computer assisted tomography [20] (for a more complete account of applications, see [21]). In seismics, the RT is known as the slant stack or the tau-P transform. The inversion procedure for the RT was rediscovered for the purpose of processing seismograms probably as early as in 1969 [9]. Heuristically, the use of the RT in signal processing can be explained as follows: a seismogram (a multidimensional signal array) can be viewed as a superposition of different events with energy concentrated along straight lines (at least locally). The RT maps these events into points thus allowing identification and separation. The RT was used by many investigators for transformation and analysis of seismograms [5]-[9], [12]-[14] and, also, for computing synthetic seismograms [ 1 11. An applicable theory of the generalized RT (integration over surfaces or curves with the presence of a weight function) was developed in [2]-[4]. Recently, inversion procedures, known in geophysical literature as migration algorithms, were cast as inversions of the causal generalized RT [22]-[24]. This makes it even more important to look carefully into the problem of exact and fast inversion of discrete versions of the RT and its generalizations. The approach we adopt in this paper can be illustrated by a parallel with the discrete Fourier transform (DFT). Similar to the DFT, the DRT is defined and studied as a transform in its own right. The DFT, of course, is used to compute direct and inverse continuous Fourier transform integrals, especially using the fast Fourier transform (FFT) algorithm. Similarly, we show how DRT can be used to compute the classical RT, and the inversion procedure for DRT can be used to invert it. Moreover. we 0096-3518/87/0200-0162$01.00